NT@UW-03-012 
BNL-NT-03/9 



Baryon Stopping and Valence Quark Distribution at Small x 



Kazunori Itakura 1 , Yuri V. Kovchegov 2 , Larry McLerran 3 and Derek Teaney 3 

m 

1 RIKEN BNL Research Center, Brookhaven National Laboratory, Bldg. 5 10 A 

Upton NY, 11973 

(N : 

Oh| 2 Department of Physics, University of Washington, Box 351560 

Seattle, WA 98195 

. 3 Nuclear Theory Group, Brookhaven National Laboratory, Bldg. 510 A 

Upton, NY 11973 

m 
> 

CO 

We argue that the amount of baryon stopping observed in the central rapidity region 
of heavy ion collisions at RHIC is proportional to the nuclear valence quark distributions 
at small x. By generalizing Mueller's dipole model to describe Reggeons we construct a 
non-linear evolution equation for the valence quark distributions at small x in the leading 
double- logarithmic approximation. The equation includes the effects of gluon saturation 
in it. The solution of the evolution equation gives a valence quark distribution function 
dn V ai/dy ~ e < - 0,4 ' °' 5 ' y . We show that this y-dependence as well as the predictions of Regge 
^ 1 theory are consistent with the net-proton rapidity distribution reported by BRAHMS. 
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I. INTRODUCTION 



One of the outstanding theoretical issues related to our understanding of hadrons is the problem of how 
distributions of quarks and gluons arise. In the last several years, remarkable progress has been made in 
our understanding of the gluon distribution at small x [1—9]. One has been able to show the existence of 
a region of a high density of gluons in a very coherent configuration, the Color Glass Condensate and to 
understand how this region joins on to the low density region [2,3,7,8,10-13]. These regions correspond to 
different values of x and Q 2 for which the gluon distribution function is measured. The region of low density 
can be understood by a combination of techniques associated with the BFKL [14] and the DGLAP [15] 
evolution equations. Within the low density region, it has been found that the geometric scaling [16] extends 
outside of the saturation region, where the correlation functions for gluons are pure powers with anomalous 
dimensions [11]. At larger Q 2 and or larger values of x, the density of gluons is even lower and results match 
on smoothly to the results of DGLAP evolution [15]. These various regions are shown in Fig. 1. (For specific 
numerical estimates see [17].) 

The essential feature of small x physics which allows for a solution for the gluon distribution is that the 
density of gluons is very large [2]. This density can be thought of as a momentum scale squared, Q 2 , where 
the subscript s refers to saturation [1,18]. Saturation here means that the density of gluons at any given 
value of x approaches a fixed limit as we go toward smaller x. For p\ < Q 2 , the gluon phase space density 
in a hadron or nucleus of radius R is [2] 
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Because the strong coupling constant is evaluated at Q 2 a ^> Aq CD , the coupling is small. The phase space 
density is large, and hence the gluons are in a condensate. The glass arises because the gluons arc described 
by classical fields which are produced by sources at higher values of x and therefore have their time scales 
Lorentz time dilated. Furthermore, the gluon distribution function is somewhat disordered in the transverse 
plane [7,8] suggesting glassy behavior. Since the coupling constant is weak, weak coupling methods may be 
combined with renormalization group techniques to solve for the properties of the Color Glass Condensate 
[4,6-8]. Because in the weak field region, the typical momentum scale is large, there is no infrared problem 
in solving BFKL or DGLAP evolution. The infrared cutoff is in fact the saturation momentum as this is the 
scale on which the gluon color charge distribution neutralizes [2,3,9,13,19]. 
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FIG. 1. The various phases of high density QCD. (For specific numerical estimates see [17]). In the Color Glass 
Condensate the gluon density is large. In the Color Quantum Fluid phase, the density is low, but correlation functions 
are power law behaved with anomalous dimensions. In the parton gas phase, the density is described by the ordinary 
evolution equations of DGLAP or BFKL equation. 



A problem currently not understood is the origin of the valence quark number distributions at small x. 
This includes both the baryon number distribution and isospin number distribution. Spin is also a valence 
effect for the distributions of quarks and gluons inside a hadron, but spin, unlike baryon number and isospin, 
can be carried by the gluons. We shall restrict our attention here to isospin and baryon number since it 
is a little simpler, although many of the techniques developed here might be applied to this case. In the 
case of baryon number, it is true that non-perturbative topological excitations of the gluon field might give 
contributions to baryon number [20-22] . For weak coupling, we expect these contributions to be small, and 
in this paper we shall not consider their effect. It would be most interesting to have a proper first principle 
computation of the magnitude of these effects within the Color Glass Condensate. An experimental analysis 
of the stopping of net electric charge versus the stopping of baryon number, which would allow one to 
disentangle the contributions of the valence quarks from the non-perturbative topological gluonic excitations 
at RHIC is under way [23]. 

We shall assume here that baryon number and isospin are carried by valence quarks [24,25]. In weak 
coupling, we should be able to compute the x and Q 2 dependence of these valence quark distribution 
functions. We will include non-perturbative aspects of the Color Glass Condensate in our computation since 
the quarks, although weakly coupled, propagate in the strong background field associated with the Color 
Glass Condensate. 

In experiment, one typically measures the valence quark distribution in the distribution of particles pro- 
duced in the final state. In this paper, we shall only consider how this is affected by the valence quark 
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distribution associated with the hadron's wavefunction. We shall not consider the effects of final state scat- 
tering. In so far as this final state scattering is due to local particle interactions, valence quantum numbers 
should spread diffusively, and cannot disperse over a wide range in x. On the other hand, novel particle 
interactions have been proposed, often involving the interpretation of baryon number as a topological exci- 
tation of gluons, which allow for transfer of baryon number [21,22]. Computation of such effects is beyond 
the scope of this paper, but they might be possible to include in a computation such as that advocated by 
Krasnitz and Venugopalan [26] and by Lappi [27]. 

In this paper, we will derive the x and px dependence of the unintegrated valence quark distribution 
functions, subject to the caveats above. To understand the essential role that the Color Glass Condensate 
plays in this computation, consider the diagrams of Fig. 2 computed by Kirschner and Lipatov [28-30] and 
by Griffiths and Ross [31]: 
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Evolution of baryon number at 
large x to small x: The double 
logarithmic approx: a ln 2 (l/x) 



FIG. 2. The ladder diagrams for the evolution of baryon number. 

These diagrams are singular, and the expansion parameter for individual diagrams is a s ln 2 (l/a;). The 
ln 2 (l/x) arises from the longitudinal and transverse phase space due to the fact that the px integral is not 
limited in the ultraviolet and is effectively cut off by the center of mass energy of the system [28-31]. 

For gluons, the expansion parameter for ladder diagrams is a s \n(l/x), and the result of summation is 
that [14] 
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For the case of valence quark distributions, we expect that the double logarithmic behavior is generated by 
l/x^ Bas , so it is not too surprising that the result of computing diagrams in Fig. 2 is [28] 
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If we take a s ~ 0.25 — 0.5, we find dn va i/dx ~ 1/a; ' 5-0 ' 7 . Such a value is typical of the Color Glass 
Condensate and the valence quark distribution, as we shall see in a later section, correctly describes the 
baryon number distribution seen at RfflC. 

The problem with summing the ladder graphs is justifying a weak coupling expansion [32]. If the coupling 
is weak, then such diagrams indeed generate the leading order contribution. On the other hand if one 
takes the ladder contribution seriously, the dominant contribution for baryon number occurs at a transverse 
momentum scale of order Aqcd, and the weak coupling methods fail. We will see that properly including 
the effects of the Color Glass Condensate generate a cutoff at px ~ Q s 3> Aqcd at asymptotically small x. 

To see how this works, consider the diagrams of Fig. 3: 
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Evolution of baryon number with - x X - 
gluon interaction in the Color 
Glass Condensate 



FIG. 3. The ladder diagrams for the evolution of baryon number including the possibility of interaction with the 
Color Glass Condensate. 

In this figure, one has a ladder sum, but the gluons can be absorbed in the Color Glass Condensate. In 
terms of traditional Feynman diagrams this implies that any number of gluonic ladder fan diagrams [1,4] 
can connect to the reggeon (quark) ladder depicted in Fig. 3. The overall interaction would look like set of 
fan diagrams consisting of a single quark ladder with any number of gluonic (BFKL [14]) ladders attached 
to it and to each other. The fan diagrams for gluons have been summed before in [1,4,6]. Here we are going 
to write down an equation summing up the fan diagrams with a single quark ladder. For small enough fey 
of the gluons in Fig. 3, the gluons do not propagate, they get absorbed in the target and the distribution in 
Pt of the quarks does not evolve. Our expectation is therefore that we get a growing distribution function 
for valence quarks only if the transverse momentum of the quarks exceed the saturation momentum. The 
coupling constant should therefore be evaluated at Q s . 

The goal of this paper is to show that this is true by constructing and solving an evolution equation for 
the valence quark distribution functions. The equation is constructed in Eq. (43) and its solution with a 
simple model for gluon evolution is given by Eq. (70). 

The resulting picture which arises for the gluon and quark distribution functions is amusing: The gluons 
have phase space density which is dominated by glue with kr < Q s and has a logarithmic divergence 
~ lnQ 3 /hr in the infrared. The valence quarks have a similar transverse phase space density, which is lower 
than the gluon one by powers of Bjorken x (see Eqs. (85), (86) and (87) and is described by the quark 
saturation scale Qi s uark = (2/3) Q s . This is shown in Fig. 4. 




FIG. 4. The phase space distribution of valence quarks and gluons as a function of pr at fixed x. 

A remarkable conclusion of this work is that for px far outside the saturation region, the valence quark 
distributions are power law behaved with non-integer powers (for fixed coupling) . This reflects the power law 
behavior with fractional anomalous dimensions found for the gluon distribution function in the geometrical 
scaling region [10,11]. It adds weight to the conjecture that there is an intermediate phase between that of 
the Color Glass Condensate and the parton gas, the Color Quantum Fluid [33]. 
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The outline of this paper is as follows: In Sect. II we construct the analog of the McLerran-Venugopalan 
model for valence quarks [2]. This illustrates with a simple model how the saturation scale regulates the 
valence quark distribution in the infrared. Next, in Sect. Ill we derive a small x evolution equation for 
the valence quark distribution within Mueller's dipole picture. This evolution equation is then solved in 
the linear regime in Sect. IV and the x-intercept is found to agree with previous results based upon the 
summation of ladder diagrams [28]. Then in Sect. V we solve the full non- linear evolution equation using a 
simple theta function model of the dipole scattering amplitude. This section therefore merges Sect. II and 
and Sect. IV illustrating how the saturation scale serves as an infrared cutoff and how quantum evolution 
changes the canonical dimensions in the McLerran-Venugopalan model. Finally, in Sect. VI we analyze the 
rapidity distribution of net-protons from the BRAHMS experiment [34] and extract a phenomenological 
intercept. This intercept is compared with the perturbative intercept of Sect. IV and the intercept expected 
from Regge theory. 



II. VALENCE QUARK DISTRIBUTION IN THE SEMI-CLASSICAL APPROXIMATION 



In this Section we construct a soft valence quark wave function of a nucleus in the quasi-classical approx- 
imation of McLerran-Venugopalan model [2] . This quasi-classical wave function will have some qualitative 
features of the full answer and will also serve as the initial condition for the evolution equation which we 
will construct below. 

Let us consider an ultrarelativistic nucleus moving in the light cone "plus" direction. Similar to [3] we 
begin by constructing soft valence quark distribution in a single nucleon at the lowest order in the coupling 
as depicted in Fig. 5. There a valence quark with momentum p splits into a quark with momentum k and 
a gluon with momentum p — k. Using the rules of light-cone perturbation theory [35] the wave function of 
the valence quark in A + = light cone gauge can be written as (see Appendix A) 

ip%x(k,p-k,z) = gT a [l + z-aX(l-z)} " (fc ^~ p ) 2 ~ ■ (4) 

where z = k + /p + , a is the quark's helicity which is conserved in the splitting since quarks are assumed to be 
massless, A and a are gluon's polarization and color and e A is the gluon's polarization vector. Transforming 
into transverse coordinate space we end up with 

^x(x 23 ,x 31 ,z) = g T a [l + z-a\{l-z)]S 2 (x 31 +zx 23 )^-^f^, (5) 

Z7T X 23 

where — x_ i — Xj with x x and x 2 the coordinates of the quark before and after the splitting and x 3 the 
transverse coordinate of the gluon (see Fig. 5). 




FIG. 5. Soft valence quark wave function of a nucleon. 



Squaring the expression in Eq. (5), summing over gluon polarizations, averaging over quark helicities and 
integrating over the initial quark position x x (in the amplitude and in the complex conjugate amplitude 
separately) we end up with 

W J 2n J Zi 2(1 -z) 2 Z-'^ AV - 23 ' n J J Zi 2n 1-z x| 3 ' K> 



5 



where we defined 



a s = 



cx s Cf 



(7) 



Zi is some initial light cone momentum fraction and Cf — N 2 C N 1 ■ Looking at the ^-dependence of Eq. (6) 
we recognize the real part of the DGLAP splitting function j qq [15]. By transforming z — > 1 — z we would 
obtain jGq splitting function as expected. 

Until now we have not imposed any restrictions on the quark's longitudinal momentum fraction z. Impos- 
ing z <C 1 for soft quarks we get x_ x = x z from Eq. (5). The soft valence quark wave function can be written 
as 



Ca(*21^) - gT a {l-aX) 



i e A • x 



21 



2tt 



(8) 



Let us define a distribution function of valence quarks in a nucleus similar to how it was done for gluons in 
[3,36]. To do that we have to multiply the wave function in Eq. (8) by its complex conjugate at a different 
transverse coordinate of the soft quark, average over helicities, sum over polarizations and integrate over 
transverse coordinate positions of the initial quark x x . The result yields 



dn val (x - y) 



dy 



LO 



1 z 
4^ 2 
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(9) 



where the brackets (...) imply averaging over nucleon's impact factors and summation over all nucleons in 
the nucleus [3,36], which brought in a factor of atomic number A. We inserted a factor of N c in Eq. (9) to 
account for N c valence quarks in a nucleon and A is some infrared cutoff of the order of Aqcd- Since we are 
interested in the number of quarks per unit rapidity we rewrote the ^-integral in Eq. (6) for small z as 

dz dz 
w dz = — z 

1 — z z 

which yielded the factor of z in Eq. (9). Eq. (9) gives us the soft valence quark distribution in a nucleus at 
the lowest order in the coupling constant. Its Fourier transform would give us unintegrated valence quark 
distribution function of a nucleus. 

Recalling that for soft gluons the corresponding lowest order wave function is [5] 



A a aX (x 21 ,z) = 2gT a — 



we obtain the LO gluon distribution function [3,36] 



dn G (x - y) 



dy 
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(x - y) 2 A 2 ' 



(10) 



(11) 



Comparing Eq. (9) to Eq. (11) we see that already at this lowest order the ratio of valence quarks to gluons 
in the small- z tail of the distribution functions is a rapidly falling function of rapidity y — In 1/z 



dn va i/dy 



dnc/dy 



■Nr. 



<~ z ~ e 



(12) 



LO 



The exact scaling with rapidity of the ratio in Eq. (12) will be modified by quantum evolution as we will see 
below. 

In the mean time let us try to construct the valence quark distribution function of a nucleus including 
the effects of multiple rescatterings [2,3,36]. Diagrammatically this is equivalent to resumming powers of 
a 2 s A 1 / z , or, equivalently, powers of x^Q 2 [3]. Following the prescription carried out for gluonic fields in [3] 
we start off in d^A^ = covariant gauge. In this gauge nucleons can not exchange gluons with each other 
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and the soft gluonic or quark fields of the nucleus are just additive with the total field of the nucleus being 
the sum of the individual fields of the nucleons. Therefore if ip^ v (x,X-) is the lowest order fermionic field 
of the valence quarks in a single nucleon in covariant gauge the field of the whole nucleus would be 

A 

^rfe*-) = x>r&*-)- ( l3 ) 

i=l 

After a gauge transformation to the A + = light cone gauge the field becomes 

A 

ipf(x,x-) = S(x,x-) ^Cfe^-) ; (14) 
i=i 

where the matrix of gauge transformation is [3,36] 

S(x,x_) = Pcxp(-igT a J d 2 b db-6(x- - b_) p a (6,6_) ln(\x-b\A)j. 
p a is a color charge density operator normalized according to 



(15) 



(p a (x,x-) p b (y,y.)) = ^- p (x,x.)S(x.-y.)S 2 (x-y)S ab (16) 
with p(x, X-) the normal nuclear density in the infinite momentum frame of the nucleus, obeying 

d 2 x dx- p(x,x-) = A. (17) 



/ 



In terms of the fermionic field ip^°(x, x_) the valence quark distribution function can be written as 

dn val (x-y) 1 z [dx-dy^ ifci ,„__„_■, / 1LCl ,1 



dy 47T 2 

= -~ J ^_^_ e ^-v-) ($?>(£,y_)S-%y-)±7 + S(x,x-)1>T(x,x-)y (18) 



4tt 2 



where now (. . .) includes averaging over longitudinal coordinates of the nucleons as well. The lowest order 
single nucleon valence quark field ip^" (x, X- ) is the same in light cone and covariant gauges and is propor- 
tional to a (5-function on the light cone, ^ v (x, X-) ~ 5(x- — Xi-) with Xj_ the light cone coordinate of the 
nucleon. Using this property of the fermion field together with Eq. (13) in Eq. (18) yields 

' / """'-~- J-|£ (^(y^S-Hy^.) l 7+ S(x, Xi .)^(M,Xi-)) • (19) 



dy 47r 2 



! = 1 



To perform the averaging in Eq. (19) one has to use the definition of S(x,X-) from Eq. (15) together with 
Eqs. (16) and (17) in the product S~ 1 (y, Xi-)S(x, Since the direction of fr_ -integration is reversed 

in S and S^ 1 we can divide all the integration into tiny slices and the first slice on the left of S~ 1 (y, Xi-) 
would include the same 6_ interval as the last slice on the right of S(x, Xi-). Since the color charge density 
correlators (16) arc local we can independently average in each slice keeping the correlators only up to 
quadratic order in p, which corresponds to the quasi-classical approximation [3]. The same procedure has 
been used for different correlators in [36,3]. Averaging of the 5*-matrices in Eq. (19) can be done independent 
of V»7+"0 since the latter term is the only one depending on the slice of the 6_-integral adjacent to (see 
the first reference in [3] for a discussion of the "last nucleon"). Performing S'-matrix averaging [36] we obtain 
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where 2L is the extent of the nucleus in the x_ -direction. The quark saturation scale is defined for a 
cylindrical nucleus as [37] 

^>Qquark2 = J ^ A ^ ^ ^ ^ 

with S± = nR 2 the cross sectional area of the nucleus with radius R and xGn(x, 1/a; 2 ) = & s ln(l/x 2 A 2 ) the 
gluon distribution in a single nucleon (see Eq. (11)). 

Noting that the correlator in Eq. (20) is, in fact, Xi- -independent and is equal to the correlator in Eq. (9), 
we can integrate over Xi- in Eq. (20). In the end we obtain the following expression for the valence quark 
structure function of a nucleus in a quasi-classical approximation 

dn va i(x-y) zN 2 S ± / -i x - v ?o* 



dy a s 7T 2 (x - y) 2 



l _ e -(*-yfQT ark V±y (22) 



Eq. (22) has an important qualitative feature in it: the valence quark distribution goes to zero at large 
transverse separations \x — y\. The effect of saturation physics on the valence quark distribution function is 
to suppress the infrared part of the distribution. This can be observed by defining the unintegrated valence 
quark distribution function 

/„„,(* = e-y,k 2 ) = -jW d 2 xe^^^. (23) 
4tt J dy 

Neglecting the x-dependence in xGn(x, 1/x 2 ) in Eq. (21) one can Fourier-transform the expression (22) 
getting 

f^k 2 ) = Z -^vU-£^\, (24) 



4 a s ir 2 \ ' Ql uark 2 

where T{y, z) is the incomplete Gamma function. The valence quark distribution (24) is much less infrared 
divergent than the Fourier transform of the lowest order expression in Eq. (9) (see Eq. (26) below). This 
property is similar to non-Abelian Weizsacker- Williams gluon distribution function [3]. One can easily see 
that in the limit of small momenta k± <C Qj uark Eq. (24) becomes 

4 a s k z k 

At very high transverse momenta k± 3> Q9 uark the valence quark distribution function (22) transformed via 
Eq. (23) maps onto f va i(x,k 2 ) given by the lowest order expression which can be obtained by either using 
Eq. (9) in Eq. (23) or by expanding Eq. (24) 



fval (•£> h. 



LO 



Z J^Aa s ^, k ± ^>Qr ark - (26) 

l k 



The quasi-classical approximation employed in this Section allowed us to derive two important features of 
the valence quark distribution at small- a;. The first feature is that, similar to the gluon distribution, multiple 
rescatterings regulate the infrared singularity as demonstrated in Eqs. (22), (24) and (25). The second feature 
is that rapidity/Bjorken x distribution of valence quarks at the classical level is just proportional to x (see 
Eq. (12)). In the following three Sections we will study how this x-dependence gets modified by inclusion of 
quantum evolution. 



III. INCLUDING NONLINEAR EVOLUTION IN THE DOUBLE LOGARITHMIC 

APPROXIMATION 

Our goal here is to describe the small-x_Bj evolution of the valence quark distribution functions including 
the effects of gluon evolution to all orders in color charge density. We consider deep inelastic scattering 
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(DIS) on a nucleus. In the dipole picture [5,38] the splitting function of a virtual photon into qq dipole is 
factorized from the subsequent evolution of the scattering cross section of the dipole on a nucleus [4] . The 
total cross section and F 2 structure function are dominated by gluon exchange at small XBj for which one 
writes 



F 2 (x Bj ,Q 2 ) = f f^^'^(x,z) 2 f d 2 bN{x,b,r 

4ir 2 a EM J 4tt J 



= In 



)• 



(27) 



Here &(x, z) describes the photon splitting into a qq dipole with transverse separation a; and moment fractions 
z and 1 — z respectively. To first order in the electromagnetic charge (see e.g. [39]) the photon splitting 
function into a qq dipole, with flavor / and electromagnetic charge (Z/e) is 



2a EM ZjN c 



[e 2 K 2 (ex) (z 2 + (1 - zf) +AQ 2 z 2 (l - zf K 2 (ex)] 



(28) 



in" t and m/ is the quark mass. Eq. (27) implicitly includes a sum over all quark 



where e 2 = z(l — z)Q 2 + m 2 
flavors. N{x,b,r) is the forward scattering amplitude at impact parameter b of the qq dipole on a nucleus 
N(x,b, t) depends upon the rapidity r = In 
approximation we substitute t = In 
such that 



A 2 

ln^- 

X B j 



where s = and z m ; n = min(z, 1 — z). In the leading log 
provided that z is not too small. N(x,b,r) is normalized 



qqA 

glue 



d 2 bN(x,b, t). 



(29) 



The quantum evolution of the gluon exchange amplitude N(x, b, t) has been resummcd to all orders in 
pomeron exchanges/color charge density [4-8] for the total cross section of a qq dipole scattering on a nucleus. 
The result was a functional differential equation [7,8] . In the large N c limit the small x Bj evolution equation 
for the forward amplitude N(x, b, r) is a nonlinear integro-differential equation [4,6] 



N(x 01 ,b,r) = 7(x 01 ,6) exp 



Aa s C F I x m 
In 



f x 01 \ a s C F f , 4a s C F ( x i . 

r H — / dr exp In (r - r ) 

V P J J 77 Jo I * VP 



X [ d2x 2^rZ2- [2^(202^+7^21^') - N{x 02 ,b+\x 2i y)N{x 12 ,b+\x w ,T% (30) 
Jp * E 02 X 12 z Z Z 

where p is some ultraviolet cutoff and 7(x 01 ,6) is the initial condition for r-evolution. In Ref. [4] , y(x 01 ,V) 
was taken in the quasi-classical Mueller-Glauber approximation [18] 

7 fe6 ) = l-e-^r^/4. (31) 

The linear part of Eq. (30) gives the BFKL equation [14,5]. Inclusion of multiple pomeron exchanges 
introduces the term quadratic in N on the right hand side of Eq. (30). In momentum space at t = and in 
the double logarithmic limit Eq. (30) reduces to GLR-MQ equation [1]. We want to construct an analogue 
of Eq. (30) for the valence quark distribution function. 

To this end, consider the difference between the structure functions of a nucleus made entirely of protons 
and a nucleus made entirely of neutrons, F 2 al = 3(F 2 — F 2 ). The gluon exchange amplitudes are identical 
for the proton and neutron and therefore F 2 al depends upon amplitudes in which a valence quark is ex- 
changed between the photon and the target. As before, we first factor the photon splitting function from the 
subsequent interaction of the qq dipole with the nucleus. Then we define the forward scattering amplitude of 
a dipole on a nucleus interacting with a single q exchange, R(x,b, z\), where z\ is the light-cone momentum 
fraction carried by the anti-quark. As in Eq. (27), the valence quark structure function can be obtained from 
R(x, 6, z\) by convoluting R with the photon splitting function 

Fr l (x Bj ,Q 2 ) = -7-¥— [ ^^<F*^%,z) 2 f d 2 b R(x,b, Zl ) . (32) 
4ir z aEM J 47r J 
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It is helpful to imagine scattering a dipole made of strange quarks (i.e. ss) on a "nucleus" made up of 
ud dipolcs in a world with three massless flavors. Then the difference between the uu and ss dipole cross 
sections is simply given by the reggeon exchange amplitude 

aZ ee aZ A *f° t A = 2jd 2 b R(x, b, Zl ) . (33) 

In the quasi-classical and double logarithmic approximations considered below only the antiquark in the 
original dipole will be responsible for the flavor exchange interaction with the nucleus. The amplitude R will 
therefore depend on the antiquark momentum fraction z\ only. 

To derive an evolution equation for R(x, b, z) we first need to construct the quasi-classical initial conditions. 
The diagrams we need to sum for that are shown in Fig. 6. For gluonic evolution (30) the initial condition 
(31) is obtained by summing up multiple rescatterings of the dipole on the nucleons in the nucleus with 
two gluons exchanged with each interacting nucleon [18]. To construct initial conditions for the valence 
quark distribution function wc modify that by replacing one of the gluonic exchanges by the quark exchange 
amplitude (see Fig. 6). Each quark exchange would bring in a suppression factor of l/s with s the center of 
mass energy. We therefore restrict ourselves to the leading term in this l/s expansion. For simplicity, we 
shall make all the quarks in the nucleus and the dipole a single flavor and therefore there are N C A valence 
quarks in the nuclear wave function which can be exchanged with the dipole. 



... i" 



FIG. 6. Forward amplitude of a qq pair interaction with the nucleus with one flavor-exchange interaction and all 
orders in gluon exchange rescatterings. 



To calculate the diagram shown in Fig. 6 we start by considering only the qq exchange part. A simple 
calculation given in Appendix B yields 

qq _ 2q s Cp f d 2 l 

val NcZiS J j2 W 

where s = (p + q) 2 is the center of mass energy per nucleon with the virtual photon carrying momentum q 
and the nucleon having momentum p. The integration over transverse momentum I of the exchanged quarks 
in Eq. (34) is logarithmically divergent both in the infrared and in the ultraviolet. We will cut it off from 
below by the same infrared cutoff A 2 that we have used in the previous Section and we will cut it off from 
above by the maximum available relevant momentum scale in the problem — the center of mass energy z\s 
of the antiquark-nucleon system. Thus Eq. (34) becomes 

- _ 2a 2 s C F 7r Zl s 

N cZl s ln A 2 ' (6b) 

The effect of multiple two gluon exchanges with different nucleons is easy to calculate and gives a factor of 
(see e.g. [39]) 

e -* 2 Qr- fc2 /4. (36) 

This factor is understood easily in the context of Glauber theory as e ~3 n<TL , where n is the nucleon density, 
a is dipole nucleon cross section, and L is the path length of the dipole in the nucleus. The factor of \ arises 
because we are calculating the amplitude while it is the amplitude squared which gives the dipole survival 
probability e~ naL . The forward amplitude R in the quasi-classical approximation is then 
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(37) 



Here a factor ( N C A) / S± has been inserted since the antiquark can be exchange with any of the (N C A) / S± 
quarks at a given impact parameter. 



FIG. 7. (A) Ladder diagram with the reggeized quarks in the f-channel and effective quark-gluon vertices considered 
in [29]; (B) Standard BFKL ladder diagram with reggeized gluons in the t-channel and effective Lipatov vertices. 

Now that we have constructed the initial condition we may start building up the small- a; evolution. Linear 
equation for a BFKL-likc ladder with quarks instead of gluons in the t-channel (see Fig. 7A) has been 
originally constructed in [29] and was subsequently studied in [30,31,40]. The resulting evolution equation 
has a peculiar property: due to the presence of quarks in the i-channel a certain part of the transverse 
momentum integration in the kernel is logarithmically divergent, similar to the integral in Eq. (34). The 
integration is effectively cut off in the ultraviolet by the center of mass energy giving an extra In s per rung 
of the ladder [29]. Therefore the ladder of Fig. 7A with quarks in the t-channel at the leading order in 
Ins effectively resums powers of the parameter a s In 2 s. This resummation is usually referred to as double 
logarithmic approximation (DLA). Note that this is in contrast to the usual BFKL ladder [14] (see Fig. 7B) 
which resums powers of a s In s and for which DLA implies resummation of a s In s In Q 2 . Of course the quark 
ladder also has an UV-safe part of the kernel resummation of which gives powers of a s Ins [30,31]. However, 
it seems a little unclear whether resummation of higher order corrections to the DLA part of the kernel would 
not give contributions parametrically of the same size as iterations of the UV-safe part of the kernel. For 
instance NLO correction to the DLA kernel would give a parametric factor of a 2 In 2 s, which is equivalent 
to double iteration of the leading logarithmic UV-safe part of the kernel (a s Ins) 2 . In what follows below we 
will work only in the double logarithmic approximation for the reggeon amplitude evolution to avoid these 
complications which may potentially require resummation of perturbation theory to all orders to find the 
leading logarithmic contribution. 

At the same time we want to resum all multiple BFKL pomeron exchanges in gluon evolution similar to 
how it was done in [4] . Each pomeron is taken in the leading logarithmic approximation giving a parametric 
contribution of the order a s p e Cc * s ln s , where p is proportional to the color charge density of the nucleus 
[2,3] and is large (p <~ a s A 1 / 3 ). Therefore, even though gluon evolution will be taken below in the leading 
logarithmic approximation while reggeon evolution will be in the DLA limit, one can see that pomeron 
exchanges of the gluon evolution are enhanced by the powers of the color charge density of the source while 
the leading logarithmic part of reggeon evolution is not. It is therefore justified to include gluon evolution in 
the leading logarithmic approximation enhanced by color charge density while keeping only the DLA part 
of the reggeon evolution. 

We want to construct an analogue of Mueller's dipole model [5] for the reggeon amplitude R. To construct 
a dipole evolution one has to take 't Hooft's large N c limit [41] in the dipole's wave function. The advantage 
of this approach is that inclusion of gluon evolution effects would become straightforward [4] . 

A single step of the evolution is the same as shown in Fig. 5. A hard valence quark splits into a soft 
quark and a hard gluon. A new color dipole 23 is created [5]. To iterate this kernel one can write down an 
evolution equation which is depicted in Fig. 8. There an initial dipole 01 either has no evolution in it at all 
(the first term on the right hand side in Fig. 8) or the antiquark 1 in the dipole splits into a gluon (double 




A 



B 
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line) and an antiquark 2 creating a new dipole 12 in addition to preexisting dipole 01. Since we are iterating 
the kernel in which the antiquark 2 is always much softer than the antiquark 1 the transverse coordinate 
of the antiquark 2 is the same on both sides of the cut. The subsequent reggeon evolution can continue 
in the dipole 12 while usual gluon evolution may happen in the dipole 01. Since in the double logarithmic 
approximation the virtual corrections are not important (since they do not give two logarithms of energy s) 
we can write down an equation using only the real part of the kernel (see e.g. [1,4]). 




FIG. 8. Evolution equation for the reggeon amplitude -R(aioi , b, Zi) in the double logarithmic approximation. As 
defined above, R(x 01 ,b, zi) is the forward amplitude of a qq dipole interacting with the nucleus by a single q exchange 
and many gluon exchanges, while N(x Q1 ,b, z\) is the forward amplitude of a qq dipole interacting with the nucleus 
by gluon exchanges only. 



Using the kernel of Eq. (6) with z = Zijz\ <C 1 we write 



a , 



, 1 min{i ia g 1 /*» 1 } dz2 d 2 X2 1 



-R(£oi>^ z i) = R o(x 01 ,b,Zi) + — / — R(x 12 ,b+ -x 2Ql z 2 ) 

Ztt J z . zi x 21 I 



x [l-Wfeoi,k*i)] (38) 

where we have switched from rapidity notation to momentum fraction z\ in the argument of N as well. (N 
depends on the rapidity of the softer quark or antiquark in the dipole [4,5], which in the case of Eq. (38) is 
the antiquark rapidity determined by Z\. ) The first term on the right hand side of Eq. (38) corresponds 
to the first term on the right hand side in Fig. 8 and represents the initial conditions given by Eq. (37). 
The transverse coordinate integral in the second term on the right hand side of Eq. (38) is logarithmically 
divergent so one has to take extra care to define the limits of the ^-integration to make it finite [29-31]. 
This property of the valence quark distribution's evolution makes it very different from the usual gluonic 
evolution. Our goal is to order dipoles in rapidity so that each newly formed dipole would have smaller 
rapidity with respect to the target than the previous dipole it was produced in. In the double logarithmic 
approximation we can assume that the antiquark in a dipole always carries a softer light cone momentum 
than the quark. Defining the rapidity of dipole 01 as Y = ln^isx^) we would get y = \n.(z 2 sx\ 2 ) for the 
rapidity of dipole 12. (Dipole 01 after splitting has roughly the same rapidity as before the splitting. Details 
of rapidity ordering discussed here for quarks are not important for the gluon evolution N .) Requiring that 
Y ^> y we get 

^2«^i4 i I ( 39 ) 

X 21 

which can also be obtained by requiring that the soft valence quark contribution dominates in the energy 
denominator. On the other hand in order to obtain logarithm of energy we need the q — > Gq splitting to 
produce a soft quark (see Eq. (6)). This translates into requirement that 

z 2 < zi- (40) 

Combining Eqs. (39) and (40) we end up with 

x 2 

2 2 <2imin{l, (41) 

X 21 
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which is the upper limit used in z 2 integral in Eq. (38). 

Eq. (38) together with Eq. (32) provide us with the way of determining the valence quark distribution 
function if the forward dipole amplitude N has been found previously from Eq. (30). Defining 

-Rfeoi^> z i) = z isR(x 01 ,b, zi) (42) 

we can rewrite Eq. (38) as 

a s /■zimin{i,x;5 1 /x! 1 } ^ g,^ l 

R(x 01 ,b, Zl ) = Ro(x 01 ,b,Z!) + 77^ — R(x 12 ,b+ -x 2Q ,z 2 ) 

2ir J z . z 2 X 21 2 



x [1-N(x {)ll b,z 1 )]. (43) 

Eq. (43) explicitly demonstrates that the z 2 integral of Eq. (43) is logarithmic and is similar in structure to 
the equation describing jet decays in the modified leading logarithmic approximation (MLLA) [42]. 

One can determine the scale for the running coupling constant in Eq. (43) following the prescription 
outlined in [43]. Since the lowest order a s N c and a s Nf running coupling corrections should come together 
to give leading order QCD beta-function, we need to calculate only the a s Nf terms to obtain the answer 
(Nf if the number of flavors) . Writing a dispersion relation for the gluon propagator in Fig. 5 one can show 
that the scale for the strong coupling constant is given by (k — zp) 2 /z, such that a s — a s [(k — zp) 2 / z] [43]. 

In the small- z limit used to obtain Eq. (43) this reduces to a s (k 2 / z). In the coordinate space for one step of 
the evolution (43) this translates into a s \z\j{z 2 x 21 )}. Therefore, running coupling effects can be included in 
Eq. (43) by replacing a s in front of the integral by a s \z\j{zi x 21 )] in the integrand. As we will see in Sect. V, 
the evolution equation (38) cuts off the infrared region with momenta less than Q s in the reggeon amplitude 
R. Similar absence of infrared diffusion was observed previously for the nonlinear evolution equation (30) in 
[44]. Therefore, together with the above argument, the scale of the coupling constant can, practically, only 
be set by either Q s or a higher momentum. In either case the coupling would be small for parametrically 
large energies. In the following, we use simply Q s for the running coupling scale. 



IV. SOLUTION OF THE LINEAR DLA EQUATION 

Let us check that the linear part of Eq. (43) is consistent with the DLA piece of the equation derived in 
[29] by explicitly solving it and comparing the result with [30,31,40]. Outside of the saturation region TV <C 1 
[10,11,4] and we can neglect it on the right hand side of Eq. (43) obtaining 



^minax^J dz2 d 2 X2 . 
2tt J z 2 x\ x 



R(x i,zi) = Mx i,zi) + 7^ J ~Z~ ~jT~ R{¥-i2i z ?) ( 44 ) 



(45) 



where fo(x 01 , z n , z-\) is obtained from Eq. (37) by putting the exponent to be equal to 1 so that 

fofeoi,«i) = a 2 a C%Trj- In ^ = r<°» (y + In — 
with 

f(°> =a^7rA (46) 

In Eq. (44) we suppressed the impact parameter dependence neglecting the shift in b on the right hand side 
which is a good approximation for a central collision of a qq pair with a large nucleus [10,11,4]. Introducing 
Laplace transform in rapidity 

R(x m , Zl ) = J ^Le^Ruten) = J ^ (z lS x 2 m Y ^fe 01 ) (47) 
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we rewrite Eq. (44) as 



where we are interested only in the azimuthally symmetric solution which is dominant at high energy. Note 
that we implicitly assume that rapidity Y > 0, so that Xq X > l/(zi s) in the first term on the right hand 
side of Eq. (44). Then ^-integration in Eq. (47) runs parallel to imaginary axis to the right of the origin. 
Defining Mellin transform 

Ruteoi) = / ^(^iA 2 )- a ^a (49) 



reduces Eq. (48) to 



which gives 



ujY 2 X(uj - A) 



f (0) (^_ A 2) 

R ^ = 9 x , x ( 51 ) 



A (wA - A 2 - Sjg.) 



Combining Eqs. (51), (49) and (47) yields 



1 duj dX , , .... 9 . 9 , x r^ ' (w 2 - A 2 ) 

= y ww^^^ (52) 

Performing the w-integration first in Eq. (52) we notice that for positive Re A the high energy asymptotics 
is dominated by the rightmost pole 

w = u,*(A) = A+!± (53) 

Eq. (52) becomes 



The integral in Eq. (54) can be done in the saddle point approximation around the saddle point A* = ^Ja s /2 
yielding 

where we switched to rapidity notation with Y = \n(zis x^). The integral in Eq. (54) can, in fact, be done 
almost exactly as will be shown in the next Section for a more general case. Using Eq. (55) in Eq. (42) 
together with Eq. (46) gives 

Q A / o \ V 4 ( v / 257-l)y JsT s 2 

i?^,^) = > 2 7r 2 ^f- -f ^(^oiA)-^ 7 6 e"^ ln ( - lA) (56) 

4 oj_ \a s J y/2a s Y 

The intercept of the reggeon in Eq. (56) is equal to 



a R = ^2h s (57) 

in agreement with [28-31,40] so that 
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R ~ e &*-VV „ e (V257-l)r_ (5g) 

To understand the transverse coordinate dependence induced by evolution let us rewrite the initial conditions 
of Eq. (37) in terms of rapidity. A simple calculation yields 

R {x,b, Zl ) = a \C%*^rxl x e- Y (y + ln-^j e -*<*r* a l*. (59) 

Comparing Eq. (59) to Eq. (56) we see that the transverse coordinate dependence of R gets modified by the 
factor of (xqi A)~^ 2 " s which also agrees with the results of [28-31,40]. We see that the small- a; evolution 
makes the valence quark distribution even more sensitive to the ultraviolet region by pushing the quarks 
toward higher transverse momenta. Overall we conclude that we have constructed a solution of the linear part 
of Eq. (38) shown in Eq. (56) and found it to be in agreement with the previous studies of the perturbative 
reggeon [28-31,40]. 



V. SOLUTION OF THE FULL NONLINEAR EQUATION WITH A SIMPLE MODEL FOR N 

Let us consider a simple model for the for the forward dipole amplitude N(x Jal ,b,T) which has correct 
qualitative features in agreement with the solution of Eq. (30) [4,6] obtained in [10,11]. Let us take 

N(x 01 ,b,r) = O^QKt) - 1) . (60) 

The saturation scale Q s changes with energy as 

Ql = K 2 e KT , (61) 

where 

r = ln^, (62) 

and k is the intercept which is usually taken to be 0.2 — 0.3 [10-12]. Eq. (60) gives us the dipole amplitude 
N which is equal to one inside the saturation region and is zero otherwise. In the spirit of the theta-function 
approximation we also model the initial conditions of Eq. (43) given in Eq. (37) by 

flo(*oi,T) = f(°»Tfl(l-4^M). (63) 

Plugging Eqs. (60) and (63) into Eq. (43) we immediately see that the solution for R(x 01 , z) is zero inside 
the saturation region and can therefore be parameterized as 

R(x 01 ,z) = 0(l-a&Q2to)£i(2oi>*)> (64) 

where we again suppress the impact parameter dependence. Substituting Eq. (64) together with Eqs. (60) 
and (63) into Eq. (43) we obtain 

MV,Y) =r* 0) p{Y + ri) + =£- dY' dr,' Ri(rf,Y') (65) 

2 Jo Jo 



where we have defined 

1 

xlQ 

and made use of the fact that 



1 = ln (66) 



Y = \n(zsx 2 ± ) = - -n (67) 
P 
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with 



P 



1 - K 



For simplicity we have also assumed that 



n n = In— = 0. 



(68) 



(69) 



Note that rj, Y > 0. 

To solve Eq. (65) we reduce it to partial differential equation, determine the boundary conditions from the 
original integral equation, and then solve the differential equation using Laplace transform methods. This 
is done in Appendix C. The complete solution is 



f(0) p 



(Y(Y + V )) 3 / 2 
(Y + T]) 3 



I 3 (2 ly /Y(Y + r,)) 



(70) 



where we have defined 7 = y . 

When Y is large and y is small, we may rewrite this formula by first employing the asymptotic expansion 
of the Bessel function and subsequently expanding as a function of y 



f(o) p 1 



=2 7 Ys 



e-T" [l + 7»7] 



7 V^W) 3 ' 2 
where we have defined the rapidity at the saturation boundary 

v 1 zs 

Ys = ln Ql- 

Re-expressing this result in terms of z and zs R(x, z) — R(x, z) = 9 (l — x^Ql) R\{x, z), we have 



(71) 



(72) 



R(x, z) 



f(0) 



7 5 / 2 y / 7TZS \Q 



z s 
2 



27 



l2\7 



l-7log(4Q«) 



log: 



3/2 



(73) 



Next we can estimate the behavior of F£ al (xBj,Q 2 ) at large Q 2 and small XBj- To determine F2 al (xBj,Q 2 ) 
we substitute Eq. (28) for photon wave function and Eq. (73) for R(z,r) into Eq. (32) for F^ 1 . For large 
Q 2 the integrand falls off very rapidly due to asymptotic descent of the modified Bessel functions. Thus 
the integral is dominated by the region when e 2 r 2 = Q 2 r 2 z(l — z) is small. This occurs when z < 

for r larger than (The region near z = 1 is suppressed relative to z — as can be seen from Eq. (73) 
.) For small values of er we approximate Ki(er) sa ^ an d neglect the corresponding contribution from 
K (er) w — log(er) to find 



Fr(x Bj ,Q<) = f-- s 1;° — 



of d 2 r friQ? dz ~ 



—Ri(r,z) . 

z 



(74) 



Next we neglect the logarithmic dependencies in the square brackets of Eq. (73) and perform the integrals 
over z and r. The result for F 2 (x, Q 2 ) is 



Fr l (x,Q 2 ) 



12\ 7 



Qt 
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(75) 



Using the leading twist relation between F$ al and the valence quark phase space distribution [45] 
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rQ 2 



we may differentiate F% al and determine the parametric behavior of the unintegrated distribution function. 

dn V ai f QsV' 1 f 1 Y 1 ^ 



dx d?k \k 2 J \x J 

A number of features warrant discussion. First, quantum evolution generates an anomalous dimension of 
—7 which enhances the valence quark distribution by a factor of ~ relative to the naive distribution 

<~ (w^J ■ Second, the intercept found for the linear case in Sect. IV , cxr = \/2a s , differs from the intercept 
found here where saturation effects were fully included, 



«« = VS- < 78 » 

The factor (1 — k)~%, reflects the fact that the saturation scale Q 2 S = A 2 e KT is serving as an infrared cutoff 
which increases with collision energy. 

VI. RELATION TO RHIC DATA ON BARYON STOPPING 

The discussion presented above addressed the issue of valence quark distributions in a single nucleus. To 
compare obtained results with the experimental data produced in heavy ion collisions at RHIC one has 
to calculate valence quark production cross section. In principle the problem can be formulated in a way 
similar to the gluon production problem in the saturation framework [46,36,26,47]. First one should solve the 
quasi-classical problem of including all multiple rescatterings in the valence quark production cross section 
[46,36,26,47] and then one should continue by including the effects of quantum evolution in the obtained 
expression [37]. The above program has been carried out for gluon production in DIS and pA collisions in 
[36,37]. However for nuclear (AA) collisions the gluon production problem complicates tremendously and still 
remains to be solved [26,47]. Here we are not going to try to solve the problem of valence quark production. 
Rather we are going to make some qualitative comparisons with the AA net-proton data assuming that the 
produced net baryons are simply proportional to the sum of the valence quarks in the incoming nuclear wave 
functions. 

The two nuclei collide with beam rapidities Yr and — Yg. At rapidity y, the valence quark content (per 
unit rapidity) of the right moving wave function is XRf va i(xR,Q 2 ), where xr = e~( YB ~ y \ Similarly the 
valence quark content of the left moving wavefunction is XRf va i(xL,Q 2 ) where xr = e~( yB+!y ) . Thus we 
expect the net baryon number to scale as 

dN net 

dk 2 dy ~ x Rfval(x R ,Q 2 s ) + X L f val {x L ,Q 2 s ) (79) 

We will parametrize, xf va i(x, Q 2 S ) with a power law ~ (x) Ar . This parametrization is motivated by Regge 
phenomenology and Ar is therefore referred to as the Reggeon intercept below. Therefore the rapidity 
dependence of the net baryons is given by 

<MjL_ „ e -A R (Y B -y) + e -A R (Y B+ y) (g0) 

dk z ay 

In the small x regime we expect the Reggeon intercept to be given by the formula 



Ah = 1-2 7 ^1-2^^. (81) 

We have fitted the observed net baryon distribution with the functional form given by Eq. (80) and 
determined the Reggeon intercept. Fig. 9 shows the fit curve with Ar = 0.47 together with the naive intercept 
Ar = 1. The naive intercept completely fails to reproduce the data. The remaining curve Ar, = 0.35 will 
be discussed shortly. 
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FIG. 9. A comparison to preliminary BRAHMS data [34] on net-proton rapidity distributions with the functional 
form given by Eq. (80). 

The value of = 0.47 should be compared to the the non-perturbative reggeon intercept. Indeed, within 
the context of Regge theory valence quark transfer at high energy is given by the 1 = 1 Regge trajectory, 
Pi fi, ^3) • • ■ • Fits to the ir~p and n + p cross sections give [48] 

<7 t V = 13.63 s 0808 + 36.02 S -° A525 (82) 
o^J = 13.63 s 0808 + 27.56 s" 04525 (83) 

The difference in these cross sections is proportional to the forward amplitude with valence quark exchange. 
Thus in Regge theory the exchange amplitude scales as, R ~ s -o.4525 ^ ^,0.4525^ xhis prediction of Regge 
theory that A# = 0.4525 seems to agree rather well with the BRAHMS preliminary net-proton rapidity 
distributions. 

We wish to compare the fitted intercept with Eq. (81). The current state of theory is not sufficient to 
compare in detail with RHIC data on net-baryon production. In addition, many of the approximations of 
the Color Glass Condensate are stretched in the kinematic window of the data. Nevertheless, it is useful to 
compare Eq. (80) to the net-baryons produced in the RHIC experiment in order to verify that this theoretical 
result is not in immediate contradiction with data. 

Employing the phenomcnological analysis of [49] , we estimate the saturation scale at the rapidity y 

Q 2 M = Q 2 s (0)e-«« (84) 

where Q 2 (0) = 2.05 GeV 2 and k « 0.25. Thus at y = 3 we find Q 2 s (y = 3) w 1.0 GeV 2 which is not a very 
large scale. Very roughly then, a s (Q s ) varies from 0.33 — 0.50 as y varies from — 3 units of rapidity and is 
not particularly small. Of course the calculations leading to Eq. (81) require a s -C 1 or Q 2 /Aq CD ^> 1. 

Further we estimate x in the kinematic window of the BRAHMS data. For the right moving nucleus we 
have roughly, xji = e~( Yn ~ v ' which varies between 0.005 — 0.1 in as y varies between — 3 units of rapidity. 
Thus xr is also not particularly small as we move forward in rapidity. 

Pressing onward, we substitute a s = 0.5 and a s = 0.33 into Eq. (81) and determine An w 0.35 and 
Ar w 0.47, respectively. Thus, provided the coupling is small, the perturbative Reggeon (Eq. (81)) also 
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reproduces the experimental intercept although with many more qualifications than the non-perturbative 
Rcggeon. 

Recent data from the dAu run from the RHIC collider probed the effects of gluon saturation in a nucleus 
[50-53]. The experiments show little or no suppression of moderate to high pt hadrons at mid-rapidity 
[50-53]. These data do not favor the prediction of high-p^ hadron suppression made in [54] based on 
extending small- a; evolution to high-p^ ~ Qs- Thus it would seem that the experiments have already ruled 
out log(l/a;) evolution in this kinematic domain. However, a number of points should be considered. First, 
baryon number evolution is governed by double logarithms a s log 2 (l/a;) as opposed to single logarithms 
a s log(l/a;) in the gluon case. Thus it is reasonable to hope that evolution effects are stronger in baryons 
than the corresponding effects in gluons. Second, multiple rescatterings [55-62] which were not considered 
in [54] introduce the Cronin enhancement [63] in R dA . The effect of small- a; evolution is to reduce this 
enhancement, eventually wiping out the Cronin effect at very high energies [64,65]. However, the effect of 
small- x evolution at moderately high energy is to somewhat reduce the Cronin maximum without eliminating 
it completely. Strictly speaking, the prediction of high-pr hadron suppression is only for px 3> Q s and 
beyond the Cronin maximum. For RHIC kinematics this means pr ^ 5.0 GcV which corresponds to rather 
large Bjorken x > 0.06. This value of x is significantly larger than the values of x relevant to net-baryon 
production at mid-rapidity x w 0.01. Therefore, it is not obvious that dAu data constrains the bulk properties 
of net-baryon production calculated here. Away from mid-rapidity (at say y = 3) x becomes larger ss 0.1, 
our calculation of net-baryon production is no longer reliable, and indeed the experiments rule out strong 
evolution effects for these values of x. Between y — and y = 3 the calculation provides a qualitative 
guide to the relative contributions of kinematic effects (leading to A# = 1) and kinematic+quantum effects 
(leading to An = 1 — yj2CFCt s /n). For these reasons we feel that the comparison with data is instructive if 
not completely justified. 



We have studied how isospin and baryon number are transported to small x. In particular we have 
studied how parton saturation affects the valence quark distribution. We first constructed the analog of 
the McLerran-Venugopalan (MV) model for valence quarks (Sect. II). The model illustrates how multiple 
rescatterings regulate the infrared singularities in the valence quark distribution. The saturation scale serves 
as an energy dependent infrared regulator as for the gluon case. For large transverse momentum, the valence 
quark phase space distribution at the quasi-classical level is 



where y = log (-). 

Employing Mueller's dipolc framework we subsequently constructed a small x evolution equation for the 
forward scattering amplitude with valence quark exchange between the dipole and the target (Sect. III). 
This equation illustrates how saturation influences the evolution of valence quark quantum numbers to small 
x. Indeed, as indicated by Eq. (38), quantum evolution stops as unitarity constraints set in. 

Next we investigated the solutions of the small x evolution equation in the linear and non-linear regions 
(Sect. IV and Sect. V). In the linear region, the solution reproduces the x intercept found previously by 
summing ladder diagrams with quark exchange [28] . Quantum evolution enhances the valence quark rapidity 
distribution relative to the MV model, changing the x dependence from x to 



Using a simple theta function model for the dipole scattering amplitude N(x,b,r), we then studied how 
parton saturation and the dynamics of the Color Glass Condensate influence the x evolution of valence quarks. 
As in the MV model, we found that the saturation scale acts as an infrared regulator which increases with 
energy as Q 2 S = A 2 e KT . The effect of quantum evolution is to change the canonical dimensions of the valence 
quark distribution. For large transverse momentum, we found 



VII. CONCLUSIONS 




(85) 



(x) 



(86) 
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dn va i ( Q 2 S \ X 7 !_ 2 / a s C F , . 

^a a ^J x ' where7 "V^TT^)- (87) 

The intercept A# = 1 — 27 is very similar to the intercept in the linear regime differing only by a factor of 
(1 — ft) -1 / 2 in 7. This difference reflects the increase of the saturation scale with energy. 

Finally we studied net-baryon rapidity distributions at RHIC and extracted a phenomcnological intercept 
from the data, A R w 0.47. This value is in line with the expectations of Regge theory [48]. For a s ~| this 
value is also in agreement with the intercept of Eq. (86) . Whether the perturbative reggeon can quantitatively 
explain the net-baryon data for phcnomcnologically relevant x remains an open question. Exciting new data 
on valence quarks at small x is coming from the d — Au run at the RHIC collider and this data will provide 
new constraints which will ultimately settle this question. 
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APPENDIX A: THE LIGHT CONE WAVE FUNCTION 

In this appendix we calculate the valence quark distribution in light cone wave function of a quark moving 
in the "plus" direction. The lowest order graph is given by Fig. 5 . Following the rules of light cone 
perturbation theory this graph is given by 

(t, _ M . mm^z . ^ _ M£ T k m . (A1) 

22fPf -T,iPi (P-k) + k -P Vk+ ^p+ 

2 

Here k + = zp + , a = ± is the helicity of the quark, pT = p° — p z = and we have divided the matrix 
element by y/p + k + as is conventional in the definition of V'cta- We are working in the light cone gauge A + =0 
where e A (fc) = (e+,e~,e A ) = (0, 2 |'f ,e A ) and e A± = Using the formulas for the matrix elements 

of spinors with definite helicities given in Table II of Ref. [35] 

^M 1+ ^M = 2 ( A2 ) 

Vp+ \J~<F 

r— 7_L r— = T 1 T j ( A ^J 

VP+ P ? 

a simple calculation reduces Eq. (Af ) to the wave function given in the text 

r*x(k,p-k,z) = gT a [l + z-o\{l-z)] ■ (A4) 



APPENDIX B: INITIAL CONDITIONS FOR QUANTUM EVOLUTION 

This appendix details the steps leading to Eq. (34) for the initial conditions of the quantum evolution. 
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FIG. 10. A forward scattering amplitude with q exchange. 

First calculate the squared amplitude shown in Fig. 10, summing over gluon hclicities and colors, and 
averaging over the quark and antiquark color and helicities 



(2iV c ) 2 ^' 1 N c f ■ 



(Bl) 



Here we have assumed Regge kinematics (see e.g. [66]) where —l 2 /s <C 1 and where we may replace I 2 with 
— Z 2 . Now integrate over the loop momentum associated with the final state phase space 

d 4 l 



I 



(2tt)< 



2Tr6 + (k()2TTd+(k 2 ) 



(B2) 



using the delta functions to eliminate the l + and I integrals. Dividing by the flux factor zp\p 2 = 2 zs we 
find the total cross section with flavor exchange 



OT - = 2C 2 a 2 f d\ 

° Val N r ZS J I 2 



(B3) 



The upper and lower limits of integral over I are discussed in the text. To relate a qq al to a cross section of a 



j2„qq-A 



dipolc on a nucleus at a given impact parameter — ^ — , we must multiply this cross section by the flux of 
valence quarks in the incoming nucleus, (N c A)/S±. Thus we find 



f]2 qq-A 

2R(x,b,z) = —^ 



S± zs 



12 



(B4) 



APPENDIX C: SOLUTION OF THE MODEL NON-LINEAR EVOLUTION EQUATION 

Our goal in this appendix is to solve the model non-linear evolution equation which replaces the scattering 
amplitude with a theta function. The model is described fully in Section V where an evolution equation for 
Ri(rj,Y) was obtained 

[Y r-ri+Y-Y' 

Mv,Y) - ? {0) p(Y + v) + -y- J a dY' j drj' Ri(r]', Y') . (CI) 

Here we have defined the following variables which are described fully in the text: r = ln^f, Q 2 = 
A 2 e KT , 77 = In ^Vi2 ) Y = ln(zsa; 2 | ), and p = j^— . n is the intercept and is usually taken to be 0.2 — 0.3. 
From these definitions we derive: ^ = In |p , and Y = ^ — 77 . To solve this integral equation we reduce it 
to a partial differential equation, deduce the boundary conditions using the original integral equation, and 
finally solve the differential equation employing Laplace transforms. 
Differentiating Eq. (CI) with respect to 77 we get 
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drj 



= f^p+^ ( Y dY'R^ + Y-Y'X) 
* Jo 



Differentiating Eq. (CI) with respect to Y we get 
dRi(v,Y) 



dY 



= r ~(0) p+ ^±P £ dr] > R l{r] > t Y) + ^f J dY'Ri(rj + Y-Y>, Y') 



Subtracting Eq. (C2) from Eq. (C3) we end up with 

dRi(v,Y) _ dRi(r], Y) & sP m 



dY 



dr] 



driR x (ri,Y). 



Differentiating Eq. (C4) with respect to 77 we get 

d 2 -Ri(r,,F) _ a 2 i?!(77,r) a sP ~ 



dY drj drj 2 
The initial conditions for Eq. (C5) are given by 

fl 1 (r ? ,y = 0) = f^p V , 

which follows from Eq. (CI), 



-fRi{r,,Y). 



dY 



= f(°) P 1 + 

y=o V 4 



a s PV 



which follows from Eq. (C3) combined with Eq. (C6), and 



dRi(v,Y) 



dY 



dRi(v,Y) 



r/=0 



dr] 



r/=0 



(C2) 



(C3) 



(C4) 



(C5) 



(C6) 



(C7) 



(C8) 



which follows from Eq. (C4). The solution of Eq. (C5) is uniquely specified by three boundary conditions 
[Eqs. (C6), (C7) and (C8)], which one can see explicitly by putting it on the lattice in 77, Y space. 
Let us search for solution of Eq. (C5) in the following form 



MV,Y) = J ^W^R X , 



(C9) 



where lo(X) is some unknown function of A and the A-integral runs parallel to the imaginary axis to the right 
of the origin. Plugging Eq. (C9) into Eq. (C5) we easily obtain 



w(A) = A + 



a s P 
2A 



so that Eq. (C9) becomes 



RiiV,Y)= /^-exp 



dX 



Rx, 



(CIO) 



(Cll) 



where R\ should be fixed by initial conditions [Eqs. (C6), (C7) and (C8)]. They translate into three 
equations, which are correspondingly 



J^,e^R x = f^prj, 



(C12) 



dX ^ ( x+ &*\ k ^ J l + ^zt 



2iri 



2A 



(C13) 
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and 



/ 



dX 



2ni 



: exp 



I« A = o. 



Eqs. (C12) and (C13) fix R x to be 



f(°)p ~ 



(C14) 



(C15) 



71=1 



where a„'s are arbitrary constants, i.e. R\ is fixed by the first two conditions up to an analytic function, 
which is specified only by the third condition [Eq. (C14)]. Next, substitute Eq. (C15) into Eq. (C14) and 
perform the contour integrals using the integral representation of the modified Bessel function 1 



I Az) = ^dr I A 1 "ex,, (A • — ) 



2ni 



(C16) 



where the contour starts at — oo above the real axis, circles origin, and returns to — oo below the real axis. 
We obtain 



= 



2r(°) 



h ( V25 s P Y 2 ) + d J h ( V25 s P Y 2 ) + a 2 h ( ^2a s p Y 2 ) + higher order I' m s. (C17) 
a s \ Z Z 



Since all of the Bessel functions in Eq. (C17) are linearly independent and depend upon Y only, the coefficient 
in front of each of I v should be 0. Enforcing this condition we end up with 



0-2 = - 



4f(0) 
paf 



for n ^ 2. 



We thus have 



Rx = 



f(Q) p 4fW 
A 2 " - ; " 



which, together with Eq. (Cll) gives 

dA 



JJi(»7,K) 



2m 



: exp 



A 2 p 2 a 2 



A' 



Integrating over A in Eq. (C20) yields the solution quoted in the text 

Y + r] 



Ri(v,Y) 



1 



V Y (Y + V) 



(C18) 



(C19) 



(C20) 



(C21) 



where we have defined 7 = J . 
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